function plotrocks14(filenumber,regionalplot)
    N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
     Nsurfaces=readdata(0,24);
     surfaces=readdata(0,25);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=1; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=32; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
Iinilcrust=find(mtype==2);
clear('mcs','mct');

oldtime=time;

for count=filenumber
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
     %P=readdata(count,23);
     Nsurfaces=readdata(count,24);
     surfaces=readdata(count,25);

    dt=time-oldtime; oldtime=time;
    x=zeros(1,S); z=zeros(1,T); 
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    figure(1); clf; cmap=jet(256); cmap(1,:)=[1 1 1]; %colormap(cmap);

    
    Iconvertedcrust=find(mtype(Iinilcrust)==1);
    mtype(Iinilcrust(Iconvertedcrust))=4;
    mtype(Ib)=mtype(Ib)+0.5;
    
    mTemp=readdata(count,15);
    [grid,xgrid,zgrid]=regbilingrid(mTemp,mx,mz,x(1),x(end),z(1),z(end),2*S,2*T);
    contvals=1200:10:1800; cont=contourc(xgrid,zgrid,grid',contvals);
    clear('mTemp');


    if regionalplot
        h1=subplot(3,1,1);
        [grid,xgrid,zgrid]=regbilingrid(mtype,mx,mz,x(1),x(end),z(1),z(end),2*S,2*T); ind0=find(grid<0.75); grid(ind0)=nan;
        imagesc(xgrid,zgrid,grid'); set(gca,'ydir','rev'); axis equal; axis([0 2100000 0 710000]); caxis([0 5.5]); colorbar; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); hold on;
        hold on;
    end

    if regionalplot h2=subplot(3,1,2); else h2=subplot(2,1,1); end
    %xmin=x(S)/2-100e3; xmax=x(S)/2+100e3; zmin=10e3; zmax=150e3;
    xmin=0; xmax=120e3; zmin=0; zmax=30e3;
    [grid,xgrid,zgrid]=regbilingrid(mtype,mx,mz,xmin,xmax,zmin,zmax,2*(xmax-xmin)/(x(end)-x(1))*S,2*(zmax-zmin)/(z(end)-z(1))*T); ind0=find(grid<0.75); grid(ind0)=nan;
    imagesc(xgrid,zgrid,grid'); set(gca,'ydir','rev'); axis equal; axis([xmin xmax zmin zmax]); caxis([0 5.5]); colorbar; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); hold on; clear('Xgrid','Zgrid','grid','mtype');
    hold on;
    
    pos=1; allypoints=cell(1,length(contvals));
    while pos<length(cont)
        step=cont(2,pos);
        i=find(contvals==cont(1,pos));
        xcont=cont(1,pos+1:pos+step); ycont=cont(2,pos+1:pos+step);
        if regionalplot plot(h1,xcont,ycont,'k','linewidth',0.1); end
        %plot(h2,xcont,ycont,'k','linewidth',0.1);
        pos=pos+step+1;
        allypoints(i)={[cell2mat(allypoints(i)) ycont]};
    end
    
    clim=[0 2];
    if regionalplot subplot(3,1,3); else subplot(2,1,2); end 
    %mstraintot=(readdata(count,18).^2+readdata(count,19).^2).^0.5;
    mstraintot=readdata(count,20);
    [grid,xgrid,zgrid]=regbilingrid(mstraintot,mx,mz,xmin,xmax,zmin,zmax,2*(xmax-xmin)/(x(end)-x(1))*S,2*(zmax-zmin)/(z(end)-z(1))*T);
    %grid=log10(grid); %grid(grid<clim(1))=clim(1);
    imagesc(xgrid,zgrid,abs(grid')); set(gca,'ydir','rev'); axis equal; caxis(clim); colorbar; axis([xmin xmax zmin zmax]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    clear('grid','mstraintot');
    print(gcf, '-djpeg100','-r200',['outfiles\rockzoomfinitestrain' num2str(count)]);
 
end
